## Tensor or Kronecker product of two matrices

from PyM import *

def hadamard_matrix(n):
    if n==0: return matrix([[1]])
    if n > 0: 
        H = hadamard_matrix(n-1)
        k = nrows(H)
        H1 = matrix(Z_,k,k)
        for i in range(k):
            for j in range(k):
                H1[i,j] = - (H[i,j])
        return stack(splice(H,H), splice(H,H1))
        
H1 = hadamard_matrix(1)
H2 = hadamard_matrix(2)
H3 = hadamard_matrix(3)


def tensor (A,B):
    (ra,ca) = shape(A)
    (rb,cb) = shape(B)
    m = ra*rb; n = ca*cb
    T = matrix(Z_,m,n)
    for i in range(m):
        for j in range(n):
            (ia,ja) = (i//rb,j//cb)
            (ib,jb) = (i%rb,j%cb)
            T[i,j] = A[ia,ja]*B[ib,jb]
    return T
            
show(tensor(H1,H1))

show(tensor(H1,H2))
show(H3)